library(MASS)
library(tidyverse)
library(Gifi)
library(readxl)
library(countrycode)
library(WDI)
library(sandwich)
library(lmtest)

rm(list=ls())

nzep <- read_excel("rawdata.xlsx") %>% 
  filter(actor_type=="Country") %>% 
  select(name, country, 
         end_target, end_target_year, end_target_status, interim_target_year, 
         reporting_mechanism, territorial_emissions, carbon_credit_offsets,
         consumption_emissions,
         gasses_coverage, 
         international_aviation, international_shipping)
        
# excluding "achieved (self-declared)" countries
nzep <- nzep %>% 
  filter(name!="Benin" & name!="Bhutan" & name!="Gabon" &
           name!="Guyana" & name!="Cambodia" & name!="Suriname") 
# excluding "no target" countries
nzep <- nzep %>% 
  filter(name!="Bolivia" & name!="Botswana" & name!="Cayman Islands" &
           name!="Georgia" & name!="Libya" & name!="Oman" &
           name!="Poland" & name!="Syrian Arab Republic") 

nzep <- nzep %>% 
  filter(name!="European Union")

colnames(nzep)

# 0) end_target
unique(nzep$end_target)
netzero_equiv <- c("Net zero", "Climate neutral", "Carbon neutral(ity)", "Zero carbon")
nzep$end_target_code <- ifelse(nzep$end_target %in% netzero_equiv, 1, 0)
unique(nzep$end_target_code)
nzep <- nzep %>% 
  filter(end_target_code==1) # excluding non-net-zero countries

write.csv(nzep, "nzep.csv")
unique(nzep$end_target)

# 1) end_target_year
unique(nzep$end_target_year)
nzep$end_target_year_code[nzep$end_target_year>2050] <- 0
nzep$end_target_year_code[nzep$end_target_year==2050] <- 1
nzep$end_target_year_code[nzep$end_target_year<2050] <- 2
unique(nzep$end_target_year_code)

# 2) end_target_status
unique(nzep$end_target_status)
nzep$end_target_status_code[is.na(nzep$end_target_status) |
                              nzep$end_target_status=="Proposed / in discussion"] <- 0
nzep$end_target_status_code[nzep$end_target_status=="Declaration / pledge"] <- 1
nzep$end_target_status_code[nzep$end_target_status=="In policy document"] <- 2
nzep$end_target_status_code[nzep$end_target_status=="In law"] <- 3
unique(nzep$end_target_status_code)

# 3) interim_target_year
unique(nzep$interim_target_year)
nzep$interim_target_year_code[is.na(nzep$interim_target_year)] <- 0
nzep$interim_target_year_code[!is.na(nzep$interim_target_year)] <- 1
unique(nzep$interim_target_year_code)

# 4) reporting_mechanism
unique(nzep$reporting_mechanism)
nzep$reporting_mechanism_code[nzep$reporting_mechanism=="No reporting mechanism"] <- 0
nzep$reporting_mechanism_code[nzep$reporting_mechanism=="Less than annual reporting"] <- 1
nzep$reporting_mechanism_code[nzep$reporting_mechanism=="Annual reporting"] <- 2
unique(nzep$reporting_mechanism_code)

# 5) territorial_emissions
unique(nzep$territorial_emissions)
nzep$territorial_emissions_code[nzep$territorial_emissions=="Not Specified" |
                                  is.na(nzep$territorial_emissions)] <- 0
nzep$territorial_emissions_code[nzep$territorial_emissions=="Yes"] <- 1
unique(nzep$territorial_emissions_code)

# 6) gasses_coverage
unique(nzep$gasses_coverage)
nzep$gasses_coverage_code[nzep$gasses_coverage=="Not Specified" |
                            nzep$gasses_coverage=="Carbon dioxide only"] <- 0
nzep$gasses_coverage_code[nzep$gasses_coverage=="Carbon dioxide and other GHGs"] <- 1
unique(nzep$gasses_coverage_code)

# 7) carbon_credit_offsets
unique(nzep$carbon_credit_offsets)
nzep$carbon_credit_offsets_code[nzep$carbon_credit_offsets=="Not Specified" |
                                  nzep$carbon_credit_offsets=="Yes" |
                                  is.na(nzep$carbon_credit_offsets)] <- 0
nzep$carbon_credit_offsets_code[nzep$carbon_credit_offsets=="No"] <- 1
unique(nzep$carbon_credit_offsets_code)

# 8) consumption_emissions
unique(nzep$consumption_emissions)
nzep$consumption_emissions_code[nzep$consumption_emissions=="Not Specified" |
                                  nzep$consumption_emissions=="No" |
                                  is.na(nzep$consumption_emissions)] <- 0
nzep$consumption_emissions_code[nzep$consumption_emissions=="Yes"] <- 1
unique(nzep$consumption_emissions_code)

# 9) international_aviation
unique(nzep$international_aviation)
nzep$international_aviation_code[nzep$international_aviation=="Not Specified" |
                                  nzep$international_aviation=="No" |
                                  is.na(nzep$international_aviation)] <- 0
nzep$international_aviation_code[nzep$international_aviation=="Yes"] <- 1
unique(nzep$international_aviation_code)

# 10) international_shipping
unique(nzep$international_shipping)
nzep$international_shipping_code[nzep$international_shipping=="Not Specified" |
                                   nzep$international_shipping=="No" |
                                   is.na(nzep$international_shipping)] <- 0
nzep$international_shipping_code[nzep$international_shipping=="Yes"] <- 1
unique(nzep$international_shipping_code)

# 1) end_target_year T
# 2) end_target_status T
# 3) interim_target_year F
# 4) reporting_mechanism T
# 5) territorial_emissions F
# 6) gasses_coverage F
# 7) carbon_credit_offsets F
# 8) consumption_emissions (ex1) F
# 9) international_aviation (ex1) F
# 10) international_shipping (ex1) F

pcdata <- nzep[,15:24] %>% as.data.frame()

fitord <- princals(pcdata, ndim=3, ordinal=c(T,T,F,T,F,F,F,F,F,F))
summary(fitord,cutoff=0.01)
nzepscore <- fitord[["objectscores"]][,1]

sum(-nzepscore<=0)

data <- data.frame(
  countryname=nzep$name,
  countrycode=nzep$country,
  nzepscore=-nzepscore
)

guess_field(data$countrycode)

### covariates

## nd-gain

ndgain <- read.csv("gain.csv")[,c(1,2,27)] # 2019
colnames(ndgain) <- c("countrycode", "countryname", "ndgain")
data <- left_join(data, ndgain, by="countrycode")

## carbon footprint

cfp <- read_excel("carbonfootprint.xlsx")
colnames(cfp) <- c("year","origin","destination","co2eq") # Ggco2eq
cfp2019 <- cfp %>% filter(year==2019)

#dPdC: domestic production, domestic consumption (self-consumption)
#dPfC: domestic production, foreign consumption (export)
#fPdC: foreign production, domestic consumption (import)

dPdC <- cfp2019[cfp2019$origin==cfp2019$destination,] %>% 
  select(origin,co2eq) %>% 
  rename(country=origin)

dPfC <- cfp2019[cfp2019$origin!=cfp2019$destination,] %>% 
  group_by(origin) %>% 
  summarize(exportedCO2eq=sum(co2eq)) %>% 
  rename(country=origin)

fPdC <- cfp2019[cfp2019$origin!=cfp2019$destination,] %>% 
  group_by(destination) %>% 
  summarize(importedCO2eq=sum(co2eq)) %>% 
  rename(country=destination)

#totalP: domestic production total, dPdC + dPfC
#exportShare: dPfC (export) / totalP

#totalC: domestic consumption total, dPdC + fPdC
#importShare: fPdC (import) / totalC

cfpdata <- left_join(dPdC, dPfC, by="country")
cfpdata <- left_join(cfpdata, fPdC, by="country")
colnames(cfpdata) <- c("country","dPdC","dPfC","fPdC")
cfpdata$totalP <- cfpdata$dPdC + cfpdata$dPfC
cfpdata$exportshare <- 100*cfpdata$dPfC / cfpdata$totalP
cfpdata$totalC <- cfpdata$dPdC + cfpdata$fPdC
cfpdata$importshare <- 100*cfpdata$fPdC / cfpdata$totalC

cfpdata <- cfpdata %>% 
  select(country, exportshare, importshare) %>% 
  rename(countryname=country)

cfpdata$countrycode <- countrycode(cfpdata$countryname, 
                                   origin="country.name", 
                                   destination="iso3c")

data <- left_join(data, cfpdata, by="countrycode")

## WDI variables
wdi <- WDI(country="all",start=2018,end=2018,
           indicator=c("EG.FEC.RNEW.ZS","EN.ATM.GHGT.KT.CE",
                       "NY.GDP.MKTP.PP.KD","NY.GDP.PCAP.PP.KD"),
           extra=T) %>% 
  select(iso3c,region,EG.FEC.RNEW.ZS,46,
         NY.GDP.MKTP.PP.KD,NY.GDP.PCAP.PP.KD)
colnames(wdi) <- c("countrycode","region","renewshare","ghg",
                   "gdp","pcgdp")

data <- left_join(data, wdi, by="countrycode")

## etrs

iea <- read.csv("iea.csv") %>% 
  filter(!is.na(treatyname)) # delete rows that do not have a treaty name

# check if all wrong NAs are corrected
sum(is.na(iea$treatyname))

# extract year from the ratification date
iea$year <- as.Date(iea$ceif3) %>% 
  substring(1,4) %>% 
  as.numeric()

# climate subject category is not coded, so let's manually code
# take treaties that have "climate" in their name:
climate <- grep("Climate", iea$treatyname)
for (i in climate){
  iea$subject_category[i] <- "Climate"
}
# check if treates are properly assigned with climate category
unique(iea$treatyname[iea$subject_category=="Climate"])

iea_country <- iea %>%  
  filter(inclusion=="MEA") %>% # only MEAs
  filter(agreement_type!="Amendment") %>% # exclude Amendment
  filter(subject_category=="Species" | subject_category=="Freshwater" |
           subject_category=="Energy" | subject_category=="Pollution" |
           subject_category=="Ocean" | subject_category=="Habitat" |
           subject_category=="Climate")%>% # exclude no categories, "General", and "Weapon" 
  filter(!is.na(year)) %>% # exclude treaties that a country did not ratify
  filter(year>1979) %>% 
  group_by(country, year) %>% summarize(n=n()) %>% mutate(count=cumsum(n)) %>% 
  group_by(country) %>% 
  summarize(iea_count=max(count))
iea_country$iea_count <- scale(iea_country$iea_count)
iea_country$countrycode <- countryname(iea_country$country, destination="iso3c")

data <- iea_country %>% 
  select(countrycode,iea_count) %>% 
  left_join(x=data, by="countrycode") %>% 
  rename(etrs=iea_count)


## civil society participation index (V-DEM, 2019, v2xcs_ccsi)

load("vdem.rdata")
csi <- vdem %>% 
  select(country_name, year, v2xcs_ccsi) %>% filter(year==2019) %>% 
  mutate(countrycode = countryname(country_name,destination="iso3c")) %>% 
  rename(csi=v2xcs_ccsi)

data <- left_join(data, csi, by="countrycode")

data$EUmember <- 0
eu_list <- c("Austria","Belgium","Bulgaria","Croatia","Cyprus","Denmark",
             "Estonia","Finland","France","Germany","Greece","Ireland",
             "Italy","Latvia","Lithuania","Luxembourg","Malta",
             "Portugal","Romania","Slovakia","Slovenia","Spain","Sweden") #23 EU countries in our sample
data$EUmember[data$countryname.x %in% eu_list] <- 1

## final data
data <- data %>% 
  select(countryname, countrycode, pca_base, pca_ex1, 
         ndgain, renewshare, importshare, exportshare, 
         ghg, gdp, pcgdp, etrs, csi, EUmember, region)
write.csv(data, "KDPdata.csv")
rawdata <- read.csv("KDPdata.csv")

##### analysis
data <- na.omit(data)
data$region <- as.factor(as.character(data$region))
data$ghg_log <- log(data$ghg+1)
data$gdp_log <- log(data$gdp+1)
data$pcgdp_log <- log(data$pcgdp+1)

model <- lm(nzepscore ~ ndgain + renewshare + importshare + exportshare +
               ghg_log + gdp_log + pcgdp_log + 
               etrs + csi + EUmember + region,
             data=data)
coeftest(model, vcov=vcov(model, type="HC4", cluster=~region))
summary(model)

## Figure 1

top10 <- data %>% slice_max(nzepscore, n=10)
bottom10 <- data %>% slice_min(nzepscore, n=10)
top10bottom10 <- rbind(top10, bottom10)


top10bottom10 %>% 
  arrange(nzepscore) %>% 
  ggplot(aes(x=reorder(countryname,nzepscore), y=nzepscore,
             fill=countryname)) +
  geom_bar(stat="identity") +
  geom_text(aes(label=round(nzepscore,3), hjust=0.5)) +
  coord_flip() +
  scale_fill_manual(values=c('Iceland'='#33D15D',
                             'Spain'='#33D15D',
                             'UK'='#33D15D',
                             
                             'Greece'='#4FE64A','Austria'='#4FE64A',
                             'Lithuania'='#4FE64A','Hungary'='#4FE64A',
                             
                             'Myanmar'='#B8E64A',
                             'France'='#CCE64A',
                             
                             'Portugal'='#E6E14A', # top 10 
                             
                             'Somalia'='#D7AB46','Bulgaria'='#D7AB46',
                             'Burkina Faso'='#D7AB46',
                             
                             'Afghanistan'='#D79846','China'='#D79846',
                             'Mozambique'='#D79846','Saudi Arabia'='#D79846',
                             
                             'Nigeria'='#D78046',
                             
                             'India'='#D76D46','Indonesia'='#D76D46',
                             
                             
                             'Bahrain'='#D74646'
  )) +
  theme_light(base_size=15) +
  ylab("Net Zero Emission Pledge Stringency") +
  xlab("Country") +
  theme(legend.position="none") +
  ggtitle("Figure 1. Top-10 and Bottom-10 countries in NZEP stringency score")


nzep %>% 
  filter(reporting_mechanism_code==1) %>% 
  nrow()

sum(data$ghg)
